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^ ■ Abstract 

^ , We introduce a new approach to evaluate the largest Lyapunov exponent of a family 

of nonnegative matrices. The method is based on using special positive homogeneous 
\^ ■ functionals on R5^, which gives iterative lower and upper bounds for the Lyapunov 

exponent. They improve previously known bounds and converge to the real value. 
The rate of converges is estimated and the efficiency of the algorithm is demonstrated 
■ on several problems from applications (in functional analysis, combinatorics, and lan- 

. guage theory) and on numerical examples with randomly generated matrices. The 

fH I method computes the Lyapunov exponent with a prescribed accuracy in relatively 

' high dimensions (up to 60). We generalize this approach to all matrices, not necessar- 

ily nonnegative, derive a new universal upper bound for the Lyapunov exponent, and 
show that such a lower bound, in general, does not exist. 

oo ; 1 Introduction 

^ . Let us consider a family A = {Ai, . . . , Am} of linear operators acting in R'^. To each opera- 

tor Aj we associate a positive number pj so that Ylf=iPj = 1- In the sequel we assume that 
O I every family of operators is equipped with a family of numbers. Consider a random product 

Xk = Ad^. ■ ■ ■ Adi, where all indices {dj} are independent and identically distributed random 
variables; each dj takes values 1, . . . , m with probabilities pi, . . . ,pm respectively. According 
, . ■ to the Furstenberg-Kesten theorem [S] the value HXfeU^/*^ converges with probability 1 to a 

. number p, which depends only on the family A, i.e., on the operators {Aj}^^ and on the 

5^ \ probabilities {pj}"!^. The number A = logp is called the largest Lyapunov exponent of the 
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family A. In this paper we do not deal with other Lyapunov exponents, and, for the sake of 
simplicity, we omit the word "largest". This number can be defined by the following limit 
formula 

A = lim - E log \\ Ad, ■ ■ ■ AA , (1) 

where E denotes the mathematical expectation. The results of this paper can be extended 
to more general matrix distributions, but we restrict ourselves to i.i.d. matrices taking values 
on a finite set A. 

We introduce a new approach for computing the Lyapunov exponent based on using 
special positive homogeneous functionals on W^. The idea is the following: for any such 
a functional / the minimal and maximal expected value of ^ log , B G A^ over all 

X e M*^, X 7^ 0, give a lower and upper bound respectively for A. For families of nonnegative 
matrices those bounds can be effectively computed and then optimized over certain families 
of functionals /, which leads to optimal bounds /3fc < A < that, under some general 
assumptions, rapidly converge to A as /c — )■ oo. For every k both au and (3^ are found by 
solving unconstrained convex minimization problems. The rate of convergence is proved to 
be at least linear in k, but in most of practical examples it is much faster. For dimensions d 
up to 50 it usually takes less than k = 12 iterations to estimate A with the relative precision 
1%. All computations take a few minutes on a standard desktop computer. 

In Section II we describe the new approach for operators with a common invariant cone, 
then in Section III we consider two special families of functionals and the corresponding 
bounds Ok and Pk- In Section IV it is shown that for nonnegative matrices both those bounds 
can be found and optimized over the corresponding families as solutions of certain convex 
minimization problems. In Theorems [2] and [3] we prove that under some general assumptions 
on matrices we have — (ik < f-, where C is an effective constant. In Section V this 
technique is extended to all matrices, without the nonnegativity condition. We derive an 
upper bound for A, which is sharper than the classical upper bound iE {log||S||2 | B e 
A^}. On the other hand. Theorem H] proved in that section shows that there are no good 
lower bounds for the Lyapunov exponent of general matrices. Finally, in Section VI we 
compute or estimate Lyapunov exponents of special families of matrices arising in problems 
of functional analysis, combinatorics, and language theory, and also report numerical results 
with randomly generated matrices. 

Lyapunov exponents of matrices have been studied in the literature in great detail due 
to many applications in probability, ergodic theory, functional analysis, combinatorics, etc. 
(see [m EHl Uni nil [71 EB US] and references therein). A special attention has been paid to 
the case of nonnegative matrices, i.e., matrices with nonnegative entries p3| [T2| 123] . The 
problem of computing or estimating the Lyapunov exponent is known to be extremely hard. 
It is even algorithmically undecidable in general [2S]- Nevertheless, there are several methods 
for approximate computation of the Lyapunov exponent that work well in most of practical 
cases. These are methods for special families of matrices arising in applications [3 EH [16] , 
for general nonnegative matrices [171 [ID], and for general families |H [5l H]. The method 
proposed in this paper for nonnegative matrices has several advantages compared with those 
previously known: 1) it produces upper and lower bounds for the Lyapunov exponent A 
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that both converge to A with a hnear rate as the number of iterations grows. So, the 
Lyapunov exponent is sandwiched between two values. The rate of convergence is estimated 
theoretically, but in practice, as we see in numerical examples, it converges much faster. 2) 
The method works equally well for high dimensions. In examples with randomly generated 
matrices of dimension c? < 50 it computes A with a relative error less than 1% within a few 
iterations. 3) We relax the assumptions on nonnegative matrices imposed in the previous 
papers on the subject. 

The most popular upper bound used in the literature is 



where is the set of all m*^ products of matrices of length k (with the corresponding 
probabilities). For any norm || • || this bound converges to A as A; — ?• oo. Usually one takes 
the Euclidean norm || • ||2. As for the lower bounds for nonnegative matrices, the most well- 
known of them comes from the results of Key [T^, which uses the same formula, but with 
an arbitrary submultiplicative functional instead of the norm. We shall see that our bounds 
are closer to the real value of A and have a guaranteed rate of convergence as — ?■ oo. The 
theoretical reasons for that are the following: 1) In our bounds, we manage to interchange 
the Expectation- and Maximum-operations, which results in a smaller upper bound and a 
larger lower bound. 2) We do not restrict ourself to an a priori fixed functional (or norm), 
but rather we show how to optimize it over a large family of functionals. 

In Section V we extend this approach to general matrices, without the nonnegativity 
assumption, and obtain an upper bound that is better than (|2]). We also prove that such a 
lower bound for general matrices does not exist. 

2 Operators with a common invariant cone: Bounds 
for the Lyapunov exponent 

Assume that all operators Ai, . . . ,Am share a common invariant cone K C W^, which is 
supposed to be convex, closed, solid, pointed, and having its apex at the origin. For any 
points x,y ^ M.'^ we write x>y if x — yEK and x — y > if x — y G int i^'. For an 
operator A we write A > if it leaves the cone K invariant. The same notation are used for 
the dual cone K* = {v eR'^ \ inf^^xiv, x) > }. 

Consider a functional f : K ^ M_)_. In the sequel we impose the following assumptions 
on /: 

1) / is positive, i.e., /(x) > 0, whenever x 7^ 0; 

2) / is homogeneous, i.e., f{tx) = tf{x) for any x E K,t E M+. 

Now we define two values -Fmin and -Fmax for any family A and for any functional /: 




(2) 



-^min(/') Al) 



inf E { log f{Ax) 



AeA] 



-fmax(/) Al) 



sup E { log f{Ax) 



AeA]. 
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We use the short notation F^i^{f,A) = -Fmin, and the same with F^ax, if the functional / 

1 

k 



and the family A are fixed. Denote also f'^^^ = r -^min(/, ^*')- Thus, 



Fj^l = I ^inl^E {log f{Bx)\Be A'], 

(k) (k) 

and similarly with F^iL- Thus, F^^^^ is the smallest expected value of the logarithm of the 
ratio f{xY^ — all X G -ft', X 7^ 0. Let us first make the following simple observation: 

Lemma 1. For every natural k and n we have 

{k + n)Fi^^^^ < kF^l + nF^l and [k + n) F^;^"^^ > k F^H + uF^^^^. 
Proof. We have 

{k + n) = sup E { log f{Bx) - log f{x) \ B e A'^'' } = 

sup E { (log f{B^B2x) - log f{B2x) ) + ( log /{B^x) - log f{x) ) \ B^ e A\ ^2 G ^" } 



sup 

xeK 



E { log fiB^B^x) - log f{B2x) \B^eA'}+E{ log /{B^x) - log f{x) \ B^ G } 



< 



sup E { log fiB.z) - log fiz) I 5i G ^'^ } + sup E { log /(^ax) - log f{x) | ^2 G } , 
which completes the proof for Fmat"'*- The proof for Fj^^^""^ is the same. □ 



Corollary 1. For an arbitrary functional f and for every k we have -Fmm ^ -^min '^^'^ 

-Tmax — -Tmax ■ 



Lemma 2. For an arbitrary functional f and for any n we have F^^ < A < -Fi^x- In 
particular, Fmin < A < Fir 



mm 

max ■ 



Proof. It suffices to prove that Fmin < A < Fmax- Then applying this inequality to the 
family A^ and taking into account that A(^") = nX{A) one obtains F^^^ < A < FmaL 

By the compactness argument it follows that there are positive constants Ci,C2 such 
that Ci||x|| < f{x) < C2||x||. Actually, these constants are respectively the minimum and 
the maximum of f{x) on the intersection of the unit sphere with the cone K. Applying 
Corollary [H we obtain for every x E K,x 



k 



jlog I 5 G ^'^1 > - Ejlog ||5x|| - log ||x|| BeA''^ > 



k 



{log^ + log/(5x)-log/(x)|5G^'=} = ;^log^+ F^^l > - log ^ + F, 
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Since ^ E | log | 5 G } — t- A and ^ log ^ + Fmin — > -^min as /c — )■ oo, we see that 
-^min ^ A. On the other hand, the same Corollary [T] implies that 



^ E I log f{Bx) \ b e A"] < F^ax 



for every x e K such that f{x) = 1. Furthermore, for each x G inti^ there is a positive 
constant C{x) such that for every operator B leaving the cone K invariant, we have > 
C(a;)||i?|| (see, for instance, [20j). Therefore, C2II-BII ■ ||a;|| > f{Bx) > CiC{x)\\B\\, and hence 



^E{log/(Sx) \ BeA'] 



-)■ A as — )■ 00 . 

Thus, A < Fmax- □ 

By Fekete's lemma for any sequence of nonnegative numbers {afej^gN such that 
{k + n) afc_|_„ < fcflfc + nan, k,n e 'N, the limit limj^^oodk exists and equals to inf^gpsiafc. 
Similarly, if {k + n) ak+n > ka^ + na„ , k,n G N, then limfc_>.oo Ofc = sup;,gj^afc. Apply- 
ing Lemma [U we see that limk^ooF^in = ^'^Pkm-^mL (denote this limit by F^J), and 
limfc_j.oo -^max = inffcgNFmax (dcuotc this limit by Fm^i). Invoking now Lemma [2] we obtain 
the following 

Proposition 1. For every family A and for every functional f we have 

p . < < A < < F C^) 

mm ^ min ^ ^ ^ max — max • \'->) 

Thus, to approximate the Lyapunov exponent one can take an arbitrary functional / and 
get the values F^i^ and -Fmax as a lower and upper bound respectively. Iterating, one obtains 
the bounds F^^^ and F^lx, which, by Corollary [H are, at least, not worse. If the two inner 
inequalities in (|3]) become equahties, then Fj^^^^ and Fmax converge from different sides to A, 
which allows us to compute A with an arbitrary prescribed accuracy. Sufficient conditions 
for that are given in Theorem [T] below. Finally, in the ideal case, when F^^^ = -Fmax? ah the 
inequalities in (|3]) become equahties. In this case we get a sharp value of A immediately, just 
by evaluating F^i^. Such "ideal" functionals / are called invariant. 

Definition 1. A functional f is called invariant for a family A if 

- f{x) + E I log f{Ax) AeA^= const W xeK\{0}. 

Thus, / is invariant if and only if F^[j^ = -Fmax- In view of Lemma [2] both these values 
equal to A. 

Corollary 2. For any invariant functional f the constant in DefinitionU\ is equal to the 
Lyapunov exponent A of A. 
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Certainly, invariant functionals do not exist for all families that have invariant cones. 
For nonnegative matrices sufficient conditions were obtained in [23] , we shall cite that result 
in Theorem A (Section IV). However, even if an invariant functional exists, it may be very 
difficult to find or to approximate. Nevertheless, as the following theorem says, the very 
existence of an invariant functional guarantees that for an arbitrary functional / the values 
Fj^^^ and Fmlx both converge to A with the linear rate. 

Theorem 1. For an arbitrary family A and for every functional f we have -FmTi = A. 

//, in addition, there is an invariant functional for this family, then for every functional f 
we have Pj^^^ = X. In this case 



-^mix ~ ^Lin ^ C k ^ , /C G M , 

where the constant C depends only on A and on f . 

Proof. There are positive constants Ci,C2 such that Ci||a;|| < f{x) < C2||a;||. Taking the 
operator norm = max||a,|j=i and using the fact that the mean of maxima is bigger 

than or equal to the maximum of means, we obtain for every x ^ K, \\x\\ = 1 

^ E I log ||5|| I S G I > r max E | log \\Bx\\ B eA''} > 
k I ' J k \\x\\=i I J 

I /-g^S, M I + ^(^-' I ^ ^ -^^ } = I I + • 

Taking limit as A; — )■ oo, we get A > FmTi- Comparing with ([3]) we see that A = FiLTi. 

Let / be an invariant functional for A. By the compactness argument, for an arbitrary 
functional f on K there are positive constants Ci, C2 such that Cif{x) < f{x) < C2f{x), x G 
K. Therefore, 

F^lif,^) > + ^ log ^ = ^ + ^l°gf- 



In the same way we show that Fmlx < A + | log and hence 

fLI - Fi-i < 2 (log C2 - log CO k-' . (4) 
Taking limit as — 00, we obtain Fj^J = F^^, which completes the proof. □ 

Thus, for every functional / we have Fmlx — )■ A. If the family ^possesses an invariant 
functional on the cone K, then for every functional / the values F^^^ and FmL converge 
from two sides to the Lyapunov exponent A, and the distance between them decays as 
C k~^. This provides a theoretical opportunity to compute the Lyapunov exponent with a 
given precision, using an arbitrary functional /. To realize this idea we need to compute 
the values Fmax and F^^j^ for large k. Each computation actually requires the resolution 
of an optimization problem, for which one needs to find a global optimum of the function 
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i^kix) = ^ E {log "^^^ I B G }, on the cone K. Therefore, the functional / should 
be chosen in a special way, to obtain the objective function ipk{x) convenient for global 
minimizing/maximizing. In the next section we define two families of functionals / (each 
depending on one d- dimensional parameter), and then, in Section IV, we apply them for the 
cone K = Wi_ (i.e., for the case of nonnegative matrices). Those functionals will allow us 
not only to evaluate the lower and upper bounds for A, but also to optimize these bounds 
over all values of the parameters. This leads to a fast algorithm for computing the Lyapunov 
exponent A of nonnegative matrices (Section IV). Even in relatively high dimensions (up 
to 60) that algorithm computes the Lyapunov exponent with a good precision (the relative 
error is less than 1%). The corresponding numerical examples from applications and some 
results with randomly generated matrices are given in Section VI. Then, making use of 
the semidefinite lifting technique, we partially extend our technique to general matrices, 
without the nonnegativity assumption. Applying a special functional / on the cone of 
positive semidefinite matrices, we obtain an upper bound for A, which is, at least, not worse 
than the usual upper bound ([2]) with the Euclidean norm (Section V). In practice it is much 
more efficient, which is confirmed by numerical examples in Section VI. As for effective lower 
bounds, it is shown in Section V that they actually do not exist for general matrices. This 
explains well-known negative theoretical results on the Lyapunov exponent computation [25] . 

Remark 1. We have seen that for every functional / the value -Fm^i actually coincides with 
the Lyapunov exponent. However, for this is, in general, not the case, unless the family 
A possesses an invariant functional. The main problem, therefore, is the lower bound for 
the Lyapunov exponent. In Section V we shall see examples of matrix families that have 

(k) 

no functionals / such that F^^j^ — )■ A as k ^ oo. That is why the existence of an invariant 
functional is crucial for deriving lower bounds that converge to the Lyapunov exponent. 



3 Two special functionals f(x) 

In this section we define two families of functionals /, which then will be applied to compute 
the Lyapunov exponent of nonnegative matrices. 

Let A be an arbitrary finite family of matrices sharing an invariant cone K. For every 
point a; > consider the functional /(■) = rrc{-) defined on K as follows: 

^x(l/) = mill |r>Oy<rx|. (5) 

Geometrically, this functional is a norm on K, whose unit ball is K (1 {x — K). If rx{y) < 1 , 
then y < X, therefore Ay < Ax for any operator A > 0, and hence r;j.{Ay) < r^lAx). 
Consequently, for this functional we have 

F^ax = max E | log ^^^^ AeA} = E (log r,{Ax) \AeA}. 
y>o { r^[y) J I I J 

Let us denote 

a{x) = E < log rx{Ax) \ A e A \ ; a = ini a{x) . 

I J X > 
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Similarly we define ak{x) = Fmlx and ak = m{x>o Applying now Lemma [21 we 

conclude that ak{x) > A for each x > and A; G N, and therefore > A. 

To obtain a lower estimate for A we take arbitrary v G int K* and consider the linear 
functional f{x) = {v,x). Again, this functional is a norm on K, whose unit ball is the 
intersection of K with a half-space. For this functional we have 



inf E 



log 



Denote 



inf E 

x&K 



log 



Similarly we define = -^mln ^^"^ l^k = inf„>o Pk ■ Lemma |2] now implies that (3k{v) < 
A for every v G intK* and k eN. 

Thus, we have the following bounds for the Lyapunov exponent of a matrix possessing 
an invariant cone: 

Pkiv) < A < (6) 

In general, they are not easy to compute. For instance, to evaluate /3(f) we need to find 
the global minimum over x G -ft' of the function 



(f , Ax) 
{v,x) 



(f , Ax) 
{v,x) 



AeA 



AeA 



/3 



sup 

v^intK' 



E 



log 



iy, Ax) 



AeA 



log (w, x) + ^ Pj log(f , A 



X] 



This function is not convex in x, it is actually quasiconcave, and hence its minimization may 
be hard. Nevertheless, we shall see in Section IV that in case K = the value Pk{v) is 
not only computable, but can be efficiently optimized over all v G inti^*, and the same is 
for ak{x). If we work with a general cone K, then it is more convenient to apply the linear 
functional f{x) = {v,x) to get not a lower bound (as (3k{v)), but the upper one. Doing so, 
we write F^ax for the functional f{x) = {v,x) and obtain 



7(f) 



max > Pi log (f , Aj x) . 

<,(v,x)=l ^ OK ^ J ; 



xSiK, {v,x) 



(7) 



The objective function iIj{x) = X^jli Pj ^'^si'^j ^jx) is concave, and hence its maximum on 
the convex set {x & K \ {v, x) = 1} can be efficiently found. The same can be done for 
every k: 



lk{v) 



max 

xGK, (v,x)=l 



E I log {v,Bx) 



Be A' 



(8) 



The shortcoming of this estimate is that it is very hard to minimize over the set v G int K*, 
even in the case K = MJ^. Nevertheless, choosing appropriate v one can obtain good upper 
bounds 7fc(f) that converge fast to A as — )■ oo. We use this bound in Section V for 
approximating Lyapunov exponents of general sets of matrices. 
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4 The Lyapunov exponent of nonnegative matrices 



We are going to see that in case K = M^, i.e., when all the operators Aj are written by 

nonnegative matrices, both estimates and /3k are efficiently computable. We only show 

here how to evaluate a and (3, since the computation of and (3^ is the same with replacing 

the family ^ of m matrices by the family A'' of all their m'' products of length k. 

We begin with a. Let us first note that rJy) = max ^. Therefore, 
^ ^^^^ i=l,...4 



a{x) — E < log ( max 

1^ V i=l,...,( 



{Ax) 

d Xi 



AeA 



Changing the variables, Xi — e"' , e R, we get 

d 



a{u) = E "I log ^ .'^^-^^ ^ 
Interchanging log and max we write 

a{u) — E < max log i 



, Uj — Ui 



AeA 



AeA 



Observe that the function log ( Ylj=i '^ij e"^ ) is convex in u. The maximum of convex 
functions is convex. Therefore, the value a is a solution of the following convex minimization 
problem: 



mf E < max log ( aj,- e 

«€Mrf I i=l,...,d ^ V ^ J 



AeA 



(9) 



Let us now compute (3. We have 

I3{v) = inf E I log {v, Ax) I A e ^ | 

x>0,(v,x)=l I I J 



A} 



on 



Thus, /9(f) is the minimal value of the concave function E | log {v, Ax) 

the simplex {x > , {v,x) = 1}. This minimal value is attained at an extreme point, i.e., 
at a vertex of the simplex. Since its vertices are the vectors of the canonical basis, we have 



/3{v) = -^^^d ( ~ + E I log {v,Aej) Ae A^^ . 

Since (f , e^) = vj and {v,Aej) = {v,a^), where is the jth column of the matrix A, we 
obtain 



(3{v) 



rnin^ ^ - log + E | log {v,a^) A e A^ ^ 



(10) 
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For any j and for any c G M the set of solutions f > of the inequality 



- log + E I log (f , a^) 
coincides with the set of solutions of the inequality 



AeA 



> 



exp 



E 



log {v, 



AeA 



l[{v,A,e,r 



> 



e V 



3 ' 



1=1 



which is convex, because the function YYi^i {"VyAiej)^' is concave in v. Hence, for any j the 

function —log f + E | log {v,a^) A e A] is quasiconcave in v, and therefore f3{v) is 

quasiconcave as well, as a minimum of quasiconcave functions. Thus, /3 is the solution of 
the following quasiconcave maximization problem: 



/3 



sup mm 

i;>0 j=^,---,d 



a 

log Vj + E I log 



AeA 



Let us remember that for each k the values ak and (3k are defined by formulas and (fTTl 
multiplied by ^ and with replacing ^ by Similarly for the values ak{x) and f3k{v). Thus 
for every vectors x,v > we have the following inequality: 



< A 



< 



keN. 



(12) 



4.1 Two conditions for nonnegative matrices 

The bounds ak{v) and (3k{x) are derived for all families of nonnegative matrices, and inequal- 
ity f[T^ always holds. The question, however, is do they provide really effective estimations 
for the Lyapunov exponent, i.e., do they converge to A as — )■ oo ? It appears that the 
answer is affirmative, whenever the matrices Aj are not "too sparse". More precisely, the 
family A satisfies the following two assumptions: 

(a) there is at least one strictly positive product of matrices from A (with repetitions 
permitted); 

(b) matrices from A do not have zero rows nor zero columns. 

Note that conditions (a) and (b) are assumed in most of papers studying random products 
of nonnegative matrices (see [26l HU [121 |l7l [23]). We shall see that condition (a) can always 
be omitted, but the situation with condition (b) is more complicated. 

Let us recall that (Theorem [T]), if there is an invariant functional for the family A, then 
for any other functional / the bounds Pj:^^^ and Fmlx converge linearly to the Lyapunov 
exponent. 

Theorem A For every family of nonnegative matrices satisfying conditions (a) and (h) 
there exists an invariant functional on . This functional is, moreover, concave and mono- 
tone on the cone K . 

Applying now Theorem [1] to the functionals /(■) = rx{-) and /(■) = (f, ■)) we arrive at 
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Theorem 2. For every family of nonnegative matrices satisfying (a) and (h), and for every 
vectors x,v > we have f3k{v) < A < ak{x) , k eN, and 

where the constant C depends on A,x and v. 

Remark 2. The constant C can be effectively estimated by entries of the matrices Aj, 
see [2lj. 

Since f3k > Pkiy) and ak < ak{x), we see that the estimates (3k and tend to A as 
well, and — < C k~^. In general, there is no need to evaluate the optimal values in 
problems (Q and ffTTj) with a good precision. To approximate A it suffices to find points x 
and V for which the difference ak{x) — Pk{v) is small, say, less than e, then the Lyapunov 
exponent A is found with the precision e. As we shall see in numerical examples, in practice 
the value — Pk decays much faster than and it is enough to take a reasonably small 
k (much smaller than 1/e) to compute the Lyapunov exponent A with the precision e. 

Remark 3. Estimates ak and Pk can be extended to any set of matrices sharing a polyhedral 
invariant cone K. If the cone K is spanned by vectors {/ij }^]^ and its dual cone K* is spanned 

by vectors {gi\f2i^ then writing formula for a{x) we replace -^^^r^ by ^-^^'^"^^ , and writing 

formula for /3(f) we replace by '"'"(^^l^^^l ■ Here it is important that both sets of vectors 

are finite, i.e., that the cone K is polyhedral. As we shall see in Theorem HI for families of 
matrices with a non-polyhedral invariant cone an effective lower bound for A may not exist 
at all. In particular, is not such a bound any more. 



4.2 Omitting condition (a) 

The lower and upper bounds Pk{v) and ak{x) respectively can be computed for every family of 
nonnegative matrices. However, to prove the convergence of these bounds to A we essentially 
used conditions (a) and (b), because Theorem A may fail without them. Moreover, there 
are simple examples showing that if at least one of the conditions (a) or (b) is not fulfilled, 
then the difference ak — Pk may not vanish as /c — )■ oo, in which case our bounds do not 
provide the Lyapunov exponent computation with a given prescribed accuracy. For every 
family A condition (b) can be, of course, checked immediately. Condition (a) looks more 
difficult to verify. Nevertheless is can be checked efficiently as well. The corresponding 
algorithm takes 2md^ arithmetic operations, where d is the dimension, and m is the number 
of matrices [22]. Thus, if both conditions (a) and (b) are satisfied, then we apply Theorem [2] 
to compute the Lyapunov exponent A. Otherwise, if at least one of them fails, one can still 
use inequality f lT^ . but now there is no guarantee that both parts converge to A as — )■ oo. 
For some families this still gives good numerical estimates for A (as for the binomial matrices 
from Subsection VI. 1 below). However, there are examples, when ak and Pk are too far from 
each other for all /c, and these bounds become useless. A question arises: is it possible to 
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obtain effective upper and lower bounds for the Lyapunov exponent (perhaps, different from 
Uk and Pk) without those two conditions? We do not know the answer for condition (b). 
Is it true that if nonnegative matrices are allowed to have zero rows and columns, then the 
Lyapunov exponent can be sandwiched between two efficiently computable bounds, whose 
difference tend to zero? 

As for condition (a), the answer is affirmative. That condition can be omitted, with a 
special modification of the upper bound afc(x). In this subsection we extend our approach 
to this case, when matrices of the family A do not necessarily have a positive product. To 
begin with, we need the following key result proved in [22]: 

Theorem B [22]. If a family of nonnegative matrices A satisfies condition (h), hut do not 

have a positive product, then one of the two following cases takes place: 

(1) A is reducible; 

(2) there is a partition of the set i7 = {1, . . . , c?} to r > 2 sets Qi, . . . , Qr, on which every 
matrix from A acts as a permutation. 

In the latter case there exists a product D of matrices from A that has a block-diagonal form: 
r strictly positive blocks corresponding to the sets Qi, . . .Qj.- 

Property (1) means that there is a nontrivial subspace of M.'^ spanned by several basis 
vectors, that is invariant for all matrices from A. Property (2) means that for every matrix 
A E A there is a permutation a of the set {1, . . . , r} such that AL^ C L^(^k^, k = 1, . . . ,r, 
where = { J^i^Uk ^i^i | ^ , z G fi^ } is a cone spanned by the vectors {cj | i G fi^}. 
To compute the Lyapunov exponent for a family not satisfying condition (a) one needs to 
consider both cases of Theorem B. 

Case 1. The family A is reducible. In this case there is a permutation of basis vectors, 
after which all matrices from A take a block upper-triangular form. This permutation can 
be found by a combinatorial algorithm that takes 0{d'^) arithmetic operations (see, for 
instance, [T5| Lemma 3.1] and references therein). Now it remains to refer to the main result 
of the work [9]: for a family of block upper-triangular matrices the Lyapunov exponent equals 
to the largest Lyapunov exponent of the blocks. Hence, in case (1) the problem of computing 
the Lyapunov exponent is reduced to several analogous problems in smaller dimensions. 

Case 2. There is a partition of the set f2, on which every matrix from A acts 
as a permutation. We first show that in this case, we can restrict our attention to the 
set L = U/^^ Li. This is a union of faces Lj of the positive orthant M^, corresponding to 
the sets fli of the partition. Clearly, AjL C L for each j = 1, . . . ,m. Consider an arbitrary 
functional / on the set L, which is positive {f{x) > 0, 7^ 0) and homogeneous 

{f{tx) = tf{x), X e L, t > 0). For this functional we define Fmin, F^ax, -^i'ln, -^i^n the 
same way as in Section II. Then we establish the following analogue of Lemma [2l 

Lemma 3. If a family A is irreducible, then for an arbitrary functional f on L and for each 
n we have Pj^^ < A < F^ax- In particular, < A < Fj^^x- 

Proof. The proof is literally the same as the proof of Lemma [2] with only one exception: in 
order to prove that -Fmax > A we need to show that there exists a vector x G L, for which 




as A; — oo . 



(13) 
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In the proof of Lemma [2] we established the existence of such a point x in the cone K, now 
we need this point in the set L. To this end we apply the main result of the work [13]: if a 
family A is irreducible and its matrices have no zero columns and rows, then every nonzero 
vector X G M.'^ satisfies (fT3i) . Thus, an arbitrary 7^ suffices. The remainder of the 

proof is the same as for Lemma [2j □ 



From theorem 4] it follows that for a family of matrices satisfying assumptions of the 
case (2) of Theorem B there exists an invariant functional / on L, for which -Fmin = -Fmax = A. 
Now, precisely as in the proof of Theorem [H we conclude that for every functional f on L 
one has 

FLt-Fj^l < Ck-\ ken, (14) 

where the constant C depends only on A and on /. To estimate the Lyapunov exponent it 
remains only to choose any convenient functional / in order to compute the values F^^^ and 

For -Fmin we again choose f{x) = {v,x) with arbitrary f > 0. It appears that for this 
functional we again have the equality Fmm = (3{v), where (3{v) is defined by ffTOl) . This is 
not obvious, because now we take minimum not over the whole set M^J_, but over a much 
narrower set L. We have 



F ■ 

J. rnir 



mm 

n=l 1 



inf E 

X^Ln , {v,x) = l 



log (f , Ax) 



AeA 



Since the minimum of a concave function E | log (f , Ax) \ A e A^ on the simplex {x e 
Ln , (f ,x) = 1} is attained at an extreme point, i.e., at a basis vector, we have 



min min E < lo, 

n=l,...,r j G Sl„ 



{v, Aej) 



[V, e,- 



AeA 



min E log 

j=l,...,d 



(f , 



AeA 



For -Fmax we again take the functional r.j.{y) = maxj=i ... ,i |^ , where a; > 0. For every n 



1, . . . , r we have max r^lAy) = max 



iAx)i 



F 



where a is our permutation. Whence, 



o-(n) 



max = max E < log max 

n=l,...,r L *en^(„) Xi 



(Ax) 



AeA 



This value will be denoted by a{x). Interchanging max and log and writing j = a{n), we 
obtain 

max E i max lor 



ax 



Xi 



AeA 



and 



afclx) 



1 f (Bx) 
— max E < max log 

k j=i,-,r L ien, Xj 



Be A' 



(15) 
(16) 



Applying now (1T4|) and Lemma [3] we obtain 
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Theorem 3. For every family of nonnegative matrices possessing property (2) from Theo- 
rem\^ and for every vectors x,v > we have 

Pkiv) < A < ak{x), keN, 

and 

where the constant C depends on A,x and v. 

Let us note that the values = sup^^Q f3k{v) and ak = i^fx>o(^k{x) can be efficiently 
found by using convex programming, because the function /3k{v) is quasiconcave in v, and 
a{x) is convex in u, where Xi = e'^\ i = 1, . . . ,d. The first assertion is proved, the proof of 
the second one is the same as for the function ak{x). Thus, both bounds from Theorem [3] can 
be optimized by parameters. In practice this significantly improves the rate of convergence 
to A. 

We see that for every family A of nonnegative matrices satisfying condition (b) there is an 
efficient method of computing the Lyapunov exponent. For families satisfying condition (a) 
the method is given by Theorem [21 for a reducible family (the case (1) of Theorem B) the 
problem is equivalent to several problems of smaller dimensions, for a family of the case (2) 
of Theorem B the method is given by Theorem [31 The optimal parameters v and x in 
the bounds can be effectively found by convex programming. This covers all families of 
nonnegative matrices without zero rows and columns. 

Corollary 3. For every family A of nonnegative matrices that have neither zero columns nor 
zero rows, and for every vector v > we have A > (3k{v) and A — /3k{v) < C , k eN. 



4.3 Is it possible to avoid condition (b) ? 

According to Theorem [1] for every family of nonnegative matrices, without any additional 
assumptions, the upper bound ak converges to the Lyapunov exponent A as — )■ cxd. For 
the lower bound /3k this is also true, provided the family A is irreducible and condition (b) 
is fulfilled. Moreover, the irreducibility assumption is not essential, because if the family is 
reducible, then the problem of the Lyapunov exponent computation is equivalent to several 
problems of smaller dimension. In turn, condition (b) is, in general, unavoidable. Indeed, 
if one of the matrices have a zero column, then the lower bound becomes trivial: (3k = — oo 
for all k. In some cases (but not always!) this difficulty can be overcome by special tricks. 
Let us describe two of them: 

1) Considering the transposed family. If all matrices of the family A have no zero columns 
(but may have zero rows), then the value Pk is different from — oo for each k, and may 
converge to A as — )■ oo, although this convergence is not theoretically guaranteed. Hence, 
if the matrices of A have zero columns, then one can compute Pk for the family A* = 
{Al, . . . , A^}. In many cases this leads to good lower bounds. 
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2) The modified estimate (3k{v). In some cases one can consider the following modified 
value for a vector v G Mf^, v ^ 0: 

l3k{v) = min ( — log + E < log {v,V) 

j=l,...,d ,Vj>0 \ I 

where is the jth column of the matrix B. In contrast to the value Pk{v), here the vector 
V may have zero components, and the minimum is taken over its nonzero components. It 
is shown easily that for every v this is a lower bound for A. If we choose v so that it has 
zeros at all positions corresponding to zero columns of matrices from A^, then Pk{v) may 
converge to A. We illustrate this trick in Section VI by a numerical example and apply it 
for estimating the Lyapunov exponent of large sparse matrices arising in one problem of the 
language theory. 

Let us stress that for both these tricks there are corresponding counterexamples, when 
they do not work. In general, if the matrices have zero columns/rows, we do not know any 
satisfactory lower bound for A. Finding such a bound can be considered as a challenging 
open problem. 



BeA']) , (17) 



5 The Lyapunov exponent of general matrices 

In this section we present an efficient upper bound for the Lyapunov exponent of general 
matrices (not necessarily nonnegative), and argue that such a lower bound does not exist. 

Let us first assume that all matrices Ai, . . . , Am share a common invariant cone K. The 
Lyapunov exponent A is defined as the limit of the value ^ E { log sup | B E A^ } 

x&K,\\x\\<l 

as k ^ oo. For every k this value is an upper bound for A. More generally, for any positive 
homogeneous functional on K the value 

^ E I log sup f{Bx) Be A''] (18) 



k 



x€KJ(x)<l 



is an upper bound for A, which by Theorem [T] converges to A as A; — oo. This upper bound 
is usually applied in the literature to estimate the Lyapunov exponent. The most popular 
functional / here is the Euclidean norm. On the other hand, the value 

Fit = \ sup E{log/(i?x) BeA^] (19) 
xeKj{x)<i ^ 

is, at least, not bigger (in most cases, actually, much smaller) than f lTS]) . because the max- 
imum of means does not exceed the mean of maxima. Therefore, the upper bound F^lx 
is closer to A than f lTSj) . However, its computation may face serious difficulties, because it 
involves finding the maximal value of the function ipkix) = E | log f{Bx) | B G A'' } 
on the set x G K,f{x) < 1. Basically, such a maximization problem can be efficiently 
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solved only in the case, when the function ipki^) is concave, or quasiconcave. Unfortu- 
nately, it does not possess this property for most of norms /(a;), including the Euclidean 
norm. That is why we use the linear functional f{x) = {v,x), for which the function 
tpkix) = E I log (w, 5a:) | B G A'^} is concave, and can be effectively maximized. Its 
maximum gives the upper bound 'jk{v) defined in ([7j) and in ([8]). As we have already men- 
tioned, the shortcoming of this upper bound is that the function 7fc(f) is not convex, and 
one is not able to efficiently minimize it over v G K*. Nevertheless, it is still possible to pick 
an arbitrary v, hoping to obtain a good bound for A. 

To extend this approach to general matrices, without the common invariant cone assump- 
tion, we apply the so-called semidefinite lifting. Let us have a family A = {Ai, . . . , Am}- To 
each matrix Aj we associate an operator Aj on the ^-j^ - dimensional space A^^ of symmetric 
d X (i-matrices defined as follows: 

AjX = AjXA*, X e Md. 

We thus obtain the family A = {Ai, . . . , Am}- All operators Aj now share a common 
invariant cone: the cone ICd of symmetric positive semidefinite d x ^-matrices. Let us recall 
that /C* = /C and that the scalar product in the space Aid is defined as {X,Y) = tr(XF). 
This is shown easily that X{A) = 2X{A)- Therefore, we can apply the upper bound 7a,. (f) 
to the family A and get the following upper bound for A(^) 

Tk{V) = ^ sup E {log ti (VBXB*) B e A' ] , (20) 

where X ^ means that the matrix X is positive semidefinite. For every positive definite 
matrix V the value Tk{V) is an upper bound for A(^), and it converges to it as A; —t- oo. 
To compute Tk{V) one needs to find the maximum of a smooth concave function on the 
intersection of the cone Kd with a hyperplane {X G Aid \ tr (VX) = 1}. This problem can be 
efficiently solved by standard tools of semidefinite programming (SDP, see, for instance, [3]). 
In many cases the most natural choice is to take V = I (the identity matrix), which yields 
the following upper bound: 

Tk{I) = ^ sup E {log tr(5X5*) Be A'}, (21) 

XyO, tr{X)=l 

Let us show that this upper bound is better that the standard upper bound with the Eu- 
clidean norm. 

Proposition 2. For every matrix family we have Tk{I) ^ E {log ||-B||2 | B G A'' } , A; G 
N. 

Proof. We have 

2kTk{I) < E (log sup tr (BXB*) B e A'' ] . (22) 

XtO, tr(X)=l ' 
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The supremum of the hnear function g{X) = tr [BXB*) is attained at an extreme point of 
the set X ^ 0, tr (X) = 1, i.e., at a rank one matrix X = YY*, where y G is a vector. 
Since tr(rr*) = \\Y\\l and tr {BYY*B*) = we see that sup^^o, tr(x)=i (5X5*) = 

supyeRd^ \\Y\\=i \\BY\\2 = \\B\\l, and whence the right hand side of inequahty fl22l) equals 
to 2{log||5||2 \ B eA''}. □ 

In Section VI we present numerical results with randomly generated matrices showing 
that the bound rfc(/) can indeed be much closer to A than the standard bound with the 
Euclidean norm. 

Now let us explain why there is no good lower bound for general families of matrices, even 
if they are non-factorable (do not have nontrivial common invariant subspaces in W^) and 
share a common invariant cone (not a polyhedral cone, when we have the lower bound (3k{v), 
see Remark E]). To avoid any trouble with the case A = —oo, we formulate the result for 
lower bounds of the Lyapunov radius p = e^. 

Theorem 4. a) There is a non-factorable pair A = {^1,^2} 0/ 2 x 2-matrices such that for 
every function ip on the set of all non- factorable pairs o/2 x 2-matrices, which is continuous 
at the point A and ^{A) < p{A) at any point A, we have ip{A) < p{A) — 1. 

b) There is a non-factorable pair A = {Ai,A2} of 3 x 3-matrices sharing an invariant 
cone K such that for every function ip on the set of all non-factorable pairs ofSx 3-matrices 
sharing the invariant cone K, which is continuous at the point A and p{A) < p{A) at any 
point A, we have p{A) < p{A) — 1. 

Proof, a) Consider a pair B = {Bi, B2}, where Bi is a rotation of the plane by the angle | 
about the origin, and B2 is the orthogonal projection onto the OX axis. It is easily shown 
that the Euclidean norm of every product of length k of matrices Bi, B2 is at least 2~^. Hence, 
p{B) > 1/2. Let now _Bi^„ be a rotation of the plane by the angle | ■ about the origin, 

and Bn = {Bi^n,B2}. Clearly, i3„ i3 as n 00. On the other hand, B2Bf'J^^B2 = 0, 
and hence p{Bn) = for each n E N. Therefore, if p>{Bn) < p{Bn) = for all n, then by 
continuity p{B) < 0. Thus^ p(i3) - (p{B) > 1/2. Now it remains to take A = {2Bi, 2B2}. 

b) It suffices to take A = {4i?i,4i?2}, where -61,-62 are the 2 x 2-matrices from the 
proof of part (a), Bi is the semidefinite lifting of the matrix Bi. Matrices of this pair is 
three-dimensional, and they share a positive definite cone /C3. □ 

Thus, there are pairs of matrices, whose Lyapunov exponent cannot be well-approximated 
by a continuous lower bound. This means that there is no algorithm, whose output contin- 
uously depends on the data, which for an arbitrary pair of matrices computes the Lyapunov 
exponent with a given precision. It should also be mentioned that, as it was shown by Blon- 
del and Tsitsiklis (1997), the problem of computing the Lyapunov exponent for matrices 
with integer entries is algorithmically undecidable [25] . 

Remark 4. The statement (b) of TheoremlUdoes not hold for 2 x 2-matrices, because every 
cone in is polyhedral (see Remark [3]). 
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Remark 5. The phenomenon that the Lyapunov exponent has good upper bounds, but 
does not have lower ones is explained by the fact that this value is an upper semicontinuous 
function on the families of matrices, but not lower semicontinuous. In the proof we used a 
pair of matrices, where the Lyapunov exponent is not lower semicontinuous. 



6 Numerical examples 

In this section, we show the efficiency of our method on several examples. We ffist analyze 
matrices drawn from an application in combinatorics, for which the exact value is known. 
Then, in Subsection VI. 2, we study an application in functional analysis. We then provide 
estimates for the Lyapunov exponent of some matrices arising in the language theory (Sub- 
section VI. 3). Finally, in Subsection VI. 4, we analyze diverse kinds of randomly generated 
matrices. As mentioned above, all optimization problems that we need to solve for comput- 
ing these estimates are convex unconstrained problems. We apply standard tools of convex 
optimization, namely, the matlab function fminunc, which uses a quasi-Newton procedure 
(the so-called BFGS-scheme), together with a cubic line search procedure. All computations 
took a few minutes on a standard desktop pc. 

As a general observation, our upper bound generally converges faster than the Euclidean 
bound (that is, the bound obtained from ([2]) with the Euclidean norm) towards the real 
value in the ffist steps of the algorithm. This allows us to have a fair upper estimate without 
having to compute too large products of matrices. On the other hand, for our lower bound, 
not only it also converges faster than the previously available lower estimates (see [E]), but 
in addition, generally these other lower bounds do not allow to reach a satisfactory accuracy 
at all. Even though they converge asymptotically towards the exact value, they appear to 
be often relatively far from it for the largest values of k that a classical computer can handle 
(say, /c = 14 or so). We graphically present our results in terms of the Lyapunov radius 
p = exp(A), because in most applications, it is the "physical" quantity that one wants to 
compute. 

6.1 The binomial triangle and generalizations 

Recently, the question of the density of ones in the nth row of Pascal's triangle has been 
studied in number theory. This question is equivalent to the number of odd coefficients in the 
polynomial (l+x)". In [7j, the authors generalize the question to the study of odd coefficients 
for the polynomial (1 -|- • ■ ■ -|- x"^)"". They show that for any m G N, there exists a set of 
matrices such that for large n, with probability one, the proportion of odd coefficients 
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(a) (b) 
Figure 1: Evolution of the different bounds for E2 (a) and Ee (b). 
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Moreover, it is shown in [7j that for these particular sets of matrices, it is possible to 
compute the Lyapunov exponent exactly. For this reason, we start this section on numerical 
examples with this application, for which it is possible to refer to the exact solution. Figure [1] 
shows how our estimates evolve in comparison to the exact value, and to the upper bound 
obtained from ([2]) with an arbitrary norm. We chose the Euclidean norm, as it most often 
performs best in practice. Observe that the family Eg does not satisfy condition (b) (the 
first matrix has three zero rows), hence the convergence of Pk towards A is not guaranteed 
theoretically. Nevertheless, both and ak converge rapidly towards the exact value of A. 
Moreover, as one can see on figure [1] (b), it may happen that the previously known lower 
bounds only give the trivial zero lower bound, while our lower bound rapidly converges 
towards the exact value. 



6.2 The regularity of de Rham curves 

An important application of the joint spectral characteristics of matrices is the exponents 
of global and local regularity of solutions of functional equations. In particular, they mea- 
sure the regularity of fractal curves, refinable functions, and wavelets. To a given refinable 
function, one can associate a set of matrices so that the Lyapunov exponent of this set is 
equal to the local regularity of the function at almost all points (in Lebesgue measure) of 
its domain. We consider the simplest example of refinable functions, the so-called de Rham 
curves, which are obtained from an arbitrary fiat polygon by successive cutting off its angles, 
when each side is divided by the same ratio w : (1 — 2u) : u, where u G (O, |) is a given 
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(a) (b) (c) 

Figure 2: Evolution of the different bounds for A1/3 (a), Ai/^ (b), and Ai/i (c) 



Our bounds 

- - Euclidean norm 
* other lower bounds 
_ our modified lower 
bound 



Figure 3: Evolution of the different bounds for the matrices corresponding to the language of 7/3-free 
words. 

parameter. The corresponding matrices of de Rham curve are given by 

- S ( ^ \ f l-2uj u 
•^"^ ~ \\cu l-2cu J ' \ cu 

For these (very small) matrices, our upper bound and the Euclidean bound both perform 
well when k is large, but for small values of k, our upper bound is much better. 

6.3 Words avoiding 7/3-powers 

Our last application comes from language Theory. It has been recently shown that the 
asymptotic growth of some specific languages can be approximated by joint spectral charac- 
teristics of matrices [161 [IS]- In |2], this idea has been applied to the so-called language of 
7/3-free words (reads as "seven thirds-free words"). These are words in which any subword 
is never repeated more than 7/3 times in some sense. See [2] for more information. The 
matrices corresponding to that language have size 227 x 227. In this reference, the authors 
analyze the joint and lower spectral radii of that set of matrices, but, in view of the large 
size of the matrices, nothing is said on the value of the Lyapunov exponent. The results 
are shown on Fig. [3l From this figure it appears that the true value lies within the interval 
[3.11,4.28]. We applied our techniques in order to estimate this value. The matrices in this 
application have zero rows and zero columns, therefore we use the modified lower bound Pk 
defined by (E 
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(a) 



(b) 



(d) 

Figure 4: Evolution (with the length of the computed products) of the different bounds for random matrices 
of dimension 5 ((a) and (b)) and 60 ((c) and (d)). In ((b) and (d)), we sparsify the matrices by independently 
putting each entry to zero with a probability 5/7; the other entries are i.i.d homogeneously between zero 
and one. 



6.4 Randomly generated matrices 
6.4.1 Nonnegative matrices 

We ran our algorithms on randomly selected nonnegative matrices. We report here the 
results of computation on two different sizes: reasonably small matrices of dimension 5, and 
large matrices of dimension 60. For both these sizes, we generated dense and sparse matrices. 
For the dense matrices, the entries are drawn uniformly and independently between zero and 
one, while for the second case, we sparsify the matrices by putting each entries to zero with 
a probability p = 5/7. For matrices of size 60, computing long products rapidly becomes 
prohibitive. As one can see on Fig. 16.41 (c) and (d), already with products of length 1 we 
have a fairly good approximation of the Lyapunov exponent, while the estimate with the 
Euclidean norm is still far from convergence. As one can see in figure (b), if the matrices 
become so sparse that zero rows can appear, then our lower bound (3k may vanish. In this 
case we apply the modified bound Pk- Table [1] summarizes the scalability of our method 
when the size of the matrices increases. For the same effort of computation, the accuracy 
seems to be somewhat independent of the dimension. 
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size 


k 


lower bound 


upper bound 


accuracy 


10 


12 


5.175 


5.205 


0.6% 


20 


12 


10.16 


10.22 


0.6% 


30 


11 


15.28 


15.35 


0.4% 


40 


11 


19.74 


19.78 


0.2% 


50 


11 


24.67 


24.83 


0.7% 



Table 1: Results of our algorithm on random pairs of matrices. For all matrices, each entry was drawn 
i.i.d. at random homogeneously between zero and one. In the second column, k indicates the length of the 
products computed in order to derive the bounds in columns 3 and 4. 



^ Our bounds ^ Our bounds 

J - - Euclidean norm 7 - ' - - Euclidean norm 




Figure 5: Evolution of the different bounds in the situation where some matrices have zero columns. The 
matrices have dimension 15, while 5 random columns are put to zero. In the right hand side, the matrices 
are randomly generated and then sparsified. 

We provide the results of computations (the lower bound f3k and the upper bound ak) 
for diverse sizes of matrices, as well as the relative accuracy obtained. We also mention the 
maximal length of the products computed in order to reach this accuracy. 

6.4.2 Nonnegative matrices with zero columns 

In Fig. we demonstrate applications of the modified lower bound f[T7l) for matrices with 
zero columns, when = — oo. 

6.4.3 Matrices with negative entries 

When the matrices do not have only nonnegative entries, no algorithm is known to squeeze 
the Lyapunov exponent between a lower and an upper bound. Of course, the upper iterative 
estimate ([2]) is always available, for any choice of the norm, but a good lower estimate 
converging to A as — )■ oo, most likely, does not exist at all (see Section V). The only 
thing that can be done in this situation is to improve the upper bound. In Proposition [2] 
we showed that for every k the upper bound Tk{I) is better than ([2]) with the Euclidean 
norm. Here we compare these bounds for randomly generated matrices, whose entries were 
all i.i.d. homogeneously between —0.5 and 0.5. We report in Table |2] the results for k = 1 
for matrices of size 10,20,30,40. In Fig. El we show the accuracy of the two upper bound 
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4 
3.5 

3 
2.5 

2 
1.5 

1 1.5 2 2.5 3 

Figure 6: Evolution of our upper bound T). (7) versus the classical Euclidean bound for a randomly generated 
30 X 30 matrix with positive and negative entries. 



size 


our upper bound 


Euclidean upper bound 


10 


1.5 


2.4 


20 


2.2 


4.2 


30 


2.7 


4.9 


40 


3 


5.7 



Table 2: Results of our algorithm compared with the Euclidean bound on random pairs 
of matrices with negative entries. For all matrices, each entry was drawn i.i.d. at random 
homogeneously between —0.5 and 0.5. The algorithms are applied directly on the sets of 
matrices (i.e. A; = 1). 




for a randomly generated 30 x 30 matrix. As one can see, is small already for k — 1, 

while for A; = 3 the Euclidean bound is still far from having converged. 
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